On the derivation of amplitude equations for nonlinear oscillators subject to arbitrary 

forcing 
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By using a generalization of the multiple scales technique we develop a method to derive amplitude 
equations for zero-dimensional forced systems. The method allows to consider either additive or 
multiplicative forcing terms and can be straightforwardly applied to the case that the forcing is white 
noise. We give examples of the use of this method to the case of the van der Pol-Duffing oscillator. 
The writing of the amplitude equations in terms of a Lyapunov potential allow us to obtain an 
analytical expression for the probability distribution function which reproduces reasonably well the 
numerical simulation results. 
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I. INTRODUCTION 

It is known that the general trends of the behavior of 
a dynamical system can be captured, in some cases of 
interest, by the so-called amplitude equations 0,01 de- 
scribing the slow dynamics of the envelope of the trajec- 
tories. One of the most interesting features of amplitude 
equations is their universality: many systems can share 
the same amplitude equation depending only on general 
symmetry considerations. Another advantage of using 
the amplitude equations is that they allow numerical in- 
tegration with a bigger integration step. A very success- 
ful technique, amongst others, for the derivation of am- 
plitude equations is that of the multiple scales method. 
This technique has been extensively applied, for instance, 
to both unforced and periodically forced nonlinear oscil- 
lators misllj]. However, for more general, e.g. non pe- 
riodic, forcing terms there is no systematic derivation of 
the amplitude equation. In this paper we develop a pos- 
sible extension of the multiple scales method to obtain 
amplitude equations for dynamical systems forced with 
general functions. As an application, we consider oscilla- 
tors which are randomly forced either by additive noise 
[H IE IE El or by multiplicative noise [E 0, 0, . 

The derivation of amplitude equations for randomly 
forced dynamical systems has been previously considered 
in the literature. In some works, dynamical systems are 
treated in a probabilistic way and noise effects on bifurca- 
tions are studied 0,0,EEEE- A different approach has 
been used in references |l7l HSj where the authors derive 
a stochastic Landau form as the amplitude equation for 
the stochastic Swift-Hohenberg model and use it to de- 
scribe the dynamics of the bifurcating solutions. A sim- 
plifying feature of this case is that the Swift-Hohenberg 
equation has only first order derivatives in time. Second 
order equations have been studied in at the level of 



probability distribution functions. Our method, being an 
straightforward extension of the multiple scales method, 
can be easily applied to dynamical systems including sec- 
ond order time derivatives. 

This paper is structured as follows: In section [n] the 
general theory is presented using the van der Pol-Duffing 
oscillator with an arbitrary forcing term as an example. 
In section ITTT1 the same oscillator is considered under the 
influence of an additive noise forcing term. In section llVl 
we present the results for the multiplicative noise case. 
Finally, the main conclusions of this paper are reviewed 
in section |V| 



II. GENERAL THEORY 

Although we believe our method to be quite general, 
for the sake of concreteness, in this section we study an 
oscillatory system, namely a van der Pol-Duffing oscilla- 
tor, with a general additive forcing term, as defined by 
the following equation for the variable x(t): 



x = e [fci(l-x 2 ) 



k 2 x 3 } 



+ */(*), (1) 



being f(t) any time dependent function and e is con- 
sidered as a small parameter. The specific cases fci = 
1, &2 = (van der Pol oscillator) and fci = 0, h-x = 1 
(Duffing oscillator) are contained in this general equa- 
tion. For f(t) = 0, the van der Pol-Duffing oscillator 
has an unstable fixed point and a stable limit cycle in 
the phase space (x, x). Equation parameters have been 
rescaled out such that the frequency of the linear oscil- 
lator, i.e. e = 0, is uiq — 1. For e > 0, the evolution 



can be written as x(t) = A{t)e 



iojat 



A(t)e 



, being 
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A(t) the slowly varying complex amplitude (A denotes 
the complex conjugate of A). 

We extend the method of multiple scales in order to be 
able to consider equation for a general function/ (t). 
Let us briefly review the method of multiple scales[l|- In 
this method, one looks for a series expansion of the time 



2 



dependent variable, x(t), of the form: 

x(t) = x (t) + exi(t) + e 2 x 2 (t) + 



(2) 



The main point is to consider different time scales: To, 
Ti, T 2 ,..., with T m — e m t as independent variables and 
hence any function of time becomes a function of the 
T '<?• 



(3) 
(4) 

(5) 
(6) 



where Dj 



a 

c)T, 



xo(t) = x (T ,Ti, . . . ) 
xi(t) = xi(T ,Ti, . . . ) 

f(t) = f(T ,Ti, ...). 
The time derivative is transformed according to 

d d 8 o d r^r^9^ 

(7) 

By substituting this expansion into l[TJl 
and equating coefficients up to order e , one obtains: 

(Dg + l)a;o - 0, (8) 
(Dg + l)a?i = fci(l-a:g)DoaJo-A!2ag-2DoI>iao + 
/(To.Ti,...)- (9) 

The solution of © is 

aoCTo.Ti, . . . ) = ^(?i)e iT ° + A(Ti) e - ab (10) 

where it has been assumed that the amplitude A(T{) of 
the sinusoidal solution depends only on the slow time 
variable T\ instead of all the sequence TijTa, Ac- 
cordingly, f(t) = /(T ,Ti). Replacing this solution in 
©, we get: 



2l„iT 



(D£ + l)xi = [i hA(l- \A\ 2 ) - 2iA' - 3k 2 A\A\ 2 ]e 
+[-ikiA(l - \A\ 2 ) + 2iA' - 3k 2 A\A\ 
-iA 3 e 3iT ° + iA 3 e~ 3lT ° + /(T , T x ) 
= g(T , Tx), (11) 

where A' denotes the derivative of A with respect to T\, 
and we define g(To,T\) as the r.h.s. of the equation. 
The amplitude equation is obtained by avoiding resonant 
terms with the frequency of the l.h.s of the equation, 

(e <ro | 5 (T 0) T 1 )) = ) (12) 

with the scalar product for the Tq variable defined as 

/>OG 

(v(T Q )\w(T )) = / v(T )w(T )dT , (13) 



The problem has been reduced to extract from f(t) = 
/(To,Ti) the resonant terms, i.e. those with a compo- 
nent of e lT ° . These terms are those giving a contribution 
different from zero in the previous equation. 

The particular case of a periodic forcing term, f(t) = 
cos(A t), with a forcing frequency A = uj?, + ere = 1 + 
ere, and a = 0(1) (soft resonant excitation) has been 
extensively studied before [l|. I n this case, the standard 
approach considers the decomposition: 



cos[(l + ae)t)} = cos(t + aet) = cos(T + aT x ) 
1 



and the result (e mT ° |e imT °) = 5 n , m . By using JTTJ and 
(|12|l one obtains 

- |A| 2 ) - 2A'] - ik 2 A\A\ 2 + (e lT "\f(T , TJ) = 0. 

(14) 



i). (15) 



When substituting this expression in Eq. H14[) for the A 
variable, the term multiplying e lT " is the only one giving 
a non-zero contribution. In other words, the contribution 
to the equation for A is the spectral component of the 
function f(t) at the frequency luq = 1. The resulting 
amplitude equation is: 

H = fex^d - l^l 2 ) + ^-A\A\ 2 - l -e^. (16) 

It is clear that the proper time scale for the variation of 
the amplitude A is given by T\ . 

For any other function f(t) the splitting of the time t 
in terms of slow, T , and fast, T\, variables might not be 
so straightforward. In general we will have to compute 
expressions of the form: 

(e iTo \e imTo f(T ,T 1 )). (17) 

Our proposal is to make the following decomposition: 

f(T , Ti) = e - (m - 1)jT V (m - 1)Tl/ 7(Ti/e), (18) 

which has the property that the contribution to the scalar 
product l|17fl is different from zero: 

(e iTo |e imTo /(T ,Ti)) = e ,(m ^ 1)Tl/e /(T x /e). (19) 

This simple substitution rule has to be applied with the 
necessary values of m as demanded in each case. No- 
tice that this decomposition gives always a non vanishing 
contribution to the scalar product. This is particularly 
interesting in the case of nearly constant functions f(t) 
for which the spectrum has a peak at u) = 0, far from 
the main peaks at ±luq. In the white noise case, our pro- 
posal is able to extract a suitable contribution from the 
flat spectrum. 

Coming back to our example, Eq. I|14|) contains a term 
of the form l|17|) with m = 0. Substitution of the ansatz 
Eq.(E3 yields: 

H = k^(l - \A\ 2 ) + l k 2 \ A \A\ 2 - i e -^f( Tl /e). 



(20) 
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This is our result for the amplitude equation of the 
van der Pol-Duffing oscillator under general forcing f(t). 
Concretely, the choice f(t) — cos[(l + ae)t] gives: 



dA 
dJ\ 



ki-(l- 
2 V 



.\A\ 2 )+ik 2 -A\Af 



-ioTi -i2T 1 /e 



(21) 



As compared to the standard result of Eq. i fTSjl . this 
equation includes an extra term that gives oscillations in 
the amplitude at the scale Ti/e. It is true that in this 
case these oscillations belong to the fast time scale To 
and could, in principle, be eliminated thus leading to the 
standard result, Eq. (|16l) . However, as shown in figure ^ 
the effect of this extra term is very small and decreases 
with increasing a. 

We now consider the forcing term f(t) = fo, constant, 
for which the standard method does not obtain any con- 
tribution. Our method immediately yields: 



ik 2 -A\A\ 2 



-e-^A/o- (22) 



One could again decide to eliminate the oscillations in 
the fast scale Ti/e. However, as shown in figure [3 the 
influence of the extra term can actually improve upon the 
predictions of the amplitude equation. In particular, it 
can describe the lack of symmetry around the zero value 
observed in the evolution of x(t). 

The above scheme can be applied to the case of mul- 
tiplicative forcing term. Again, we consider a specific 
example: the van der Pol-Duffing equation with a linear 
multiplicative term of the form 



+ X = e [fci(l - x 2 ) x - k 2 x 3 } + ex f(t). 



(23) 



The multiple scales ansatz @ leads in the first order 
approximation to the form (|1(J|I . The next order yields 



D 2 qXi 



h(l - x Q )D x 



Xi 

2\ 



= e lT ° 



dA 



-2i^-f + ikxA - ikxA\A\ 2 

Oil 



k 2 xl-2D D 1 xo+x f(To,T 1 ) 
3k 2 \A\ 2 A + 



/(T ,Ti)A] 
r) A 

e~ lT ° [2i— - ik x A + ikiA\A\ 2 - 3k 2 \A\ 2 A - 
oT\ 

f(T Q ,Ti)A] 



„3iT< 



°A 6 (iki + k 2 ) + e- illa A i {ik 



gi{T ,Ti). 



(24) 



The nullity of the scalar product (e lT °\gi(T Q , Ti)) = 
involves two contributions of the form l|17|l : one with 
m = 1 and another with m — — 1. These contributions 
are computed according to the general rule (JT5J. The 
amplitude equation obtained is 

— = ki-A(l-\A\ 2 ) +l k 2 6 -A\A\ 2 - -Af(Ti/e) 



(25) 



Let us be more specific and consider f(t) = cos(2t). 
The amplitude equation is: 

H = ki 1 -A(l~\A\ 2 ) + i k 2 3 -A\A\ 2 

-il-ilcos^-iie-^A. (26) 

Again the last two terms, belonging to the temporal scale 
T , could be discarded reobtaining the amplitude equa- 
tion derived after replacing directly cos(2£) = cos(2To) in 
Eq. (J23: 

^ = ki l -A(l - |A| 2 ) + i k 2 ^A\A\ 2 ~ % -A. (27) 

In figure|3|we compare the results of Eqs. (|23|l . I|26(l and 
l|27(l . Observe that both amplitude equations faithfully 
follow the envelope of the x(t) variable. 



III. ADDITIVE NOISE 

We apply the above developed technique to the case 
of an additive noise forcing term, while in a following 
section we will consider multiplicative noise. We take 
equation with a real noise term f(t) = £(i) with zero 
mean value and temporal correlations (£(£)£(£')) = 5(t — 
t'). The general amplitude equation (|2T1|l adopts now the 
form: 



V^(Tl 



(28) 

with (£(Ti)£(T{)) = 6{Ti - T{), and it has been used 
that £(Ti/e) = ^/et;(Ti). The statistical properties that 
follow from this equation are contained in the correspond- 
ing Fokker-Planck equation for the time evolution of the 
probability density function P(A,t) 20] . It is possible 
to obtain an approximate analytical expression in the 
steady state P s t(A). We first write down the real and 
imaginary parts of the amplitude equation 

= ^AAl-(A 2 r + A 2 )}-k 2 3 -MA 2 r + A 2 ) 



1 



MTi/e)V~e £(Ti), 



(29) 



dAi ki 
dJ\ ~ T 



Ai[l - (A 2 r + A 2 )} + k 2 -A r {Al + A 2 ) 



1 



cos{Ti/e)V~e t{Ti), 



(30) 



with A = A r +iAi. It can be shown that the deterministic 
terms are that of a non-relaxational potential flow with 
the Lyapunov function plT |22| : 



V(A) - ki 



\A\ 2 , \A\* 



(31) 
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Using this fact, one can obtain an approximate ana- 
lytical expression for the stationary probability distri- 
bution P st (A). To this end, we simplify the stochas- 
tic set of equations H3U|) by substituting the cosine 
and sine functions by its root mean square value 
(sin 2 ^)) 1 / 2 = (cos 2 ^)) 1 / 2 = 1/V2 and considering that 
(sin(i) cos(i)) = 0, effectively introducing different noise 
terms for the real and imaginary parts. This yields: 

where r\ is a complex noise term with zero mean value and 
correlations (r}(Ti) fj(T{)) = 8{T X -T[). It turns out that 
this new noise terms satisfy the fluctuation-dissipation 
relational], and the stationary probability distribution 
is given in terms of the Lyapunov potential as 

Pst(A) cx exp(—V(A)/e), (33) 

where e = e/16. 

We have simulated the dynamics of Eq. (JJJ with 
f(t) = £(t) following standard stohcastic integration 
methods as given in |2lJ . The results are compared with 
the corresponding amplitude equation, i.e. (|30|l . in Fig. 
0J We observe that the amplitude equation does not fit 
the actual of x(t) for a given realization. 

However, by taking many realizations, the averaged am- 
plitude does fit the maxima of the average values. As an 
evidence, we show in Fig. [5] the histogram for the x ma x 
is compared to the histogram for the amplitude. In solid 
line, it appears the approximate theoretical expression 
(|33|l . The concordance between solid line and triangles 
reveals the validity of the approximation which allowed 
us to derive |(32J| from l|5U)). 



IV. MULTIPLICATIVE NOISE 

In this section we consider the van der Pol oscillator 
with a multiplicative noise term of the form f(t) = x£(t). 
The corresponding amplitude equation l|25[l becomes: 

H = k 1 1 -A(l-\A\*) + l k 2 3 -A\Af 

-Ue^m) - l -Ae- 2iT ^V^(Tx). (34) 

In this case of multiplicative noise terms, the time traces 
are similar to the ones observed for additive noise terms, 
see Fig. EJ 

In order to be able to compute the stationary proba- 
bility distribution, we now introduce modulus and phase 



variables: A — Re^, the resulting equations are 
dR 1 , 9 , 

-^(T 1 / f ) S in(2[r 1 / £ + ffi, (35) 
§- x =^-\^h) (36) 
™eCTi/e) cosppye + fl). (37) 

Again, we approximate the sine and cosine terms by its 
root mean square value. In this way, the equation for the 
modulus does not contain the phase variable: 

This e quat ion can be analyzed using the general methods 
|2lll23l . l24l ] and it can be shown that, despite the fact that 
the Lyapunov function does not satisfy the fluctuation 
dissipation relation, Ea. l|33(l still holds asymptotically in 
the limit e — > giving: 

(Ah \ 
--i[-i? 2 + ^-]J . (39) 

This result is compared in Fig. {7\ with the numerical 
simulations. 

We would like to end this section by commenting that 
equation 138|l has the same structure that Eq. (35) of 
reference [jjj, in the sense that the noise term £(Ti) ap- 
pears multiplying the variable R. In that reference, the 
result is obtained by working at the level of the probabil- 
ity distribution function. Our method, which we believe 
to be simpler and more straightforward, works directly 
at the level of the (stochastic) dynamical equation. 

V. CONCLUSIONS 

In this work, we develop a simple, general method to 
obtain amplitude equations for a large variety of dynam- 
ical systems forced by an arbitrary time dependent func- 
tion, fit), including noises, both additive and multiplica- 
tive. Our method is based on the multiple scales analysis 
combined with a recipe to extract the resonant terms for 
arbitrary forcing functions. As a representative example, 
we analyze in detail the van der Pol-Duffing oscillator. 

In the deterministic case, our method is able to re- 
produce standard results. It only differs from them in 
the presence of some extra terms in the amplitude equa- 
tions. We show that those terms are either very small or 
improve upon the predictions of the equation. For exam- 
ple, it can incorporate the asymmetries observed in the 
simulations, as in Fig. 2. 

In the stochastic additive case, f(t) = £(t), our method 
is able to satisfactorily capture its contribution to the am- 
plitude equation. A simple set of approximations allow 
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us to obtain analytically the stationary probability distri- 
bution P st (A) in terms of a Lyapunov potential function. 

Similar conclusions can be drawn for the multiplica- 
tive noise, f(t) = x£(t). Again, our method allows us to 
obtain easily the stochastic amplitude equation. An ap- 
proximation, valid now in the limit of small noise, yields 
the stationary distribution. 
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FIG. 1: Time evolution of x(t) in a van der Pol oscillator Q 
(oscillating thin solid line) forced with f(t) — cos(l + at). In 
order to extract the envelope from the complex amplitude A 
we also plot ±27? where A — Re 1 ^ in two cases: the solid line 
is the result of method 12 U . the dashed line comes from the 
standard method llbH . Notice that both methods give almost 
indistinguishable results. The values of the parameters are 
e — 0.1, fci = 1, &2 = 0. The initial conditions: x(0) = 1, 
x(0) = are used henceforth in all the examples given. The 
different plots are: a) a = 0; b) a = 1; c) a = 10; d) a = 500. 
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FIG. 2: Time evolution of x(t) in a van der Pol oscillator 
(oscillating thin solid line) forced with f(t) = 2.0. We 
also plot the envelope ±2R(t) as coming from Eq. (solid 
line). The parameters are e = 0.1, ki = 1, fa = 0. Notice 
that the fit of the envelope to the trajectory of x(t) worsens if 
we neglect the last term in Eq. 1221 , as shown by the dashed 
line, specially in the negative values of x(t). 
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FIG. 3: Time evolution of x(t) in a van der Pol-Duffing 
oscillator 1231 (oscillating thin solid line) with multiplicative 
forcing xf(t) = xcos(2t), together with the envelope ±2R(t) 
coming both from 1261 and 1271 (they are indistinguishable 
at the scale of the figure). The values of the parameters are 
fci = 1, fa = 1, e = 0.01. 




FIG. 4: Time evolution of x(t) in a van der Pol-Duffing 
oscillator Jy (oscillating thin solid line) forced with additive 
white noise fit) = and the envelope ±2R(t) coming from 
1301 (solid line) and from the approximated expression 1321 
(dashed line). The parameters are ki = 1, fa = 1. e = 
0.1. Notice that the individual trajectories are not exactly 
approximated by the amplitude equations, while the mean 
values would be (see next figure). 
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FIG. 5: The diamonds correspond to the histogram for the 
maxima of x for a van der Pol-Duffing oscillator in the same 
case of additive white noise forcing and parameters than in 
the previous figure 0] The triangles are the histogram of the 
envelope 2R(t) as obtained from the amplitude equation 1301 . 
The solid line is the probability distribution function as given 
by the approximate analytical result of Eq. I|330 ■ 
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FIG. 6: Time evolution x(i) for a van der Pol-Duffing os- 
cillator (thin solid oscillating line) with a multiplicative noise 
term xf(t) = x£(t), and the envelope ±2R(t) coming from 
the amplitude equation 1341 (solid line). The parameters are 
fei = 1, k 2 = 1, e = 0.1. 




FIG. 7: The diamonds correspond to the histogram for the 
maxima of x for a van der Pol-Duffing oscillator in the same 
case of multiplicative white noise forcing and parameters than 
in the previous figure The triangles are the histogram of 
the envelope 2R(t) as obtained with the amplitude equation 
l|34|l . The solid line is the probability distribution function as 
given by the approximate analytical result of Eq. 139H . 



